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ABSTRACT 

A method of efficiently computing turbulent compressible flow over complex two-dimensional 
configurations is presented. The method makes use of fully unstructured meshes throughout 
the entire flow-field, thus enabling the treatment of arbitrarily complex geometries and the use 
of adaptive meshing techniques throughout both viscous and inviscid regions of the flow-field. 
Mesh generation is based on a locally mapped Delaunay technique in order to generate 
unstructured meshes with highly-stretched elements in the viscous regions. The flow equations 
are discretized using a finite-element Navier-Stokes solver, and rapid convergence to steady- 
state is achieved using an unstructured multigrid algorithm. Turbulence modeling is performed 
using an inexpensive algebraic model, implemented for use on unstructured and adaptive 
meshes. Compressible turbulent flow solutions about multiple-element airfoil geometries are 
computed and compared with experimental data. 
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ence and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 

Although unstructured mesh techniques have been employed extensively in fields such as 
solid modeling and computational structural mechanics for many years, their use in the field of 
computational fluid dynamics (CFD) constitutes a relatively recent phenomenon. This situation 
is probably due to the large overheads generally incurred with unstructured mesh techniques, 
which can become unacceptable when coupled with the large computational requirements of 
many CFD problems. The advantages of unstructured meshes lie in the ability they afford for 
flexibly discretizing arbitrarily complex geometries, and in the ease with which they lend them- 
selves to adaptive meshing techniques, which can be employed to accurately resolve complex 
flows in an efficient manner. In recent years, much progress has been made in developing 
more sophisticated unstructured mesh generation strategies for computational fluid dynamics 
problems [1,2,3] as well as in the development of novel and efficient flow-solution algorithms 
[4,5,6, 7], However, much of this effort has been directed at two- and three-dimensional invis- 
cid flow problems. The solution of high-Reynolds-number viscous flows, which are of much 
greater practical interest, introduces additional complications with regards to mesh generation 
and turbulence modeling. Most attempts at solving viscous flows using unstructured meshes 
have resorted to hybrid structured-unstructured meshes, where a thin structured mesh is placed 
in the boundary-layer and wake regions, and an unstructured mesh is employed in the regions 
of inviscid flow [8,9]. This type of compromise, however, limits the flexibility of the original 
unstructured approach. Geometries with close tolerances, where confluent wakes and boundary 
layers occur, may prove difficult to discretize with a hybrid approach, and the task of pierform- 
ing adaptive meshing throughout the viscous and inviscid regions of flow can be considerably 
more complex. 

This papier describes a method for computing compressible turbulent viscous flows about 
arbitrary two-dimensional configurations using fully unstructured meshes and incorporating 
adaptive meshing techniques. The generation of meshes with highly-stretched triangular ele- 
ments in the boundary-layer and wake regions is accomplished with a method based on a 
modified Delaunay triangulation technique [10]. The full Navier-Stokes equations are discre- 
tized and solved for on these meshes using an efficient finite-element solver which converges 
rapidly to steady-state using an unstructured multigrid strategy [11]. Turbulence modeling is 
achieved using an inexpensive algebraic model which has been devised specifically for use on 
unstructured and adaptive meshes [12]. 

2. MESH GENERATION 
2.1. Initial Mesh Generation 

The initial unstructured mesh is generated in three essentially independent stages. First, a 
distribution of mesh pioints and associated stretching vectors are generated throughout the flow 
field. These pioints are then joined together in a manner influenced by the local stretching 
values to form a set on non-overlapping triangular elements which completely fill the domain. 
The resulting mesh is then smoothed by slightly repositioning the mesh pioints according to an 
elliptic smoothing operator discretized on the new mesh. 

Although adaptive meshing techniques can be relied upon to increase the mesh resolution 
in regions of strong flow gradients, a good initial mesh-point distribution is essential in order 
to initially capture all the salient features of the flow, and to reduce the number of flow- 
solution adaptivity cycles required to converge the accuracy of the solution process. This is 
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especially true in the case of high-Reynolds-number viscous flows where a highly stretched, 
densely packed distribution is required to efficiently resolve the thin shear layers. The initial 
distribution of mesh points is generated in a geometry-adaptive manner. Points are positioned 
along all flow-field boundaries and along wake lines in a manner which is governed by the 
local rate of change of slope of the boundary curve, and the thickness (height) of the elements 
to be generated at the boundary. The effect of this distribution method is to cluster points in 
regions of high boundary curvature and sharp comers, where large flow gradients are expected, 
as well as to reduce the stretching of the mesh in such regions. 

A distribution of interior points is then constructed by generating a series of local hyper- 
bolic structured meshes for each boundary curve or wake line using the previously generated 
boundary-point distribution as an initial condition, and prescribing a normal spacing of points 
as a function of the Reynolds number of the flow to be computed [13]. The union of all the 
points contained in the various overlapping hyperbolic meshes is then used as the basis for the 
triangulation. A value of stretching is required at each mesh point, and this may be computed 
from the local ratio of the streamwise to normal point spacings in the local structured meshes, 
as shown in Figure 1. Initial point clustering in wake regions is achieved by drawing fictitious 
boundaries which approximate the estimated positions of the wakes. The position of these 
approximate wake lines is obtained by solving the corresponding inviscid flow problem using 
an inexpensive panel-method solution. t 

A point-filtering operation, based on the aspect ratios of the structured local hyperbolic 
mesh cells, is also employed to remove excessive points in the regions of inviscid flow. The 
structured-mesh cell aspect-ratios, as measured by the ratio of streamwise to normal mesh spac- 
ing, are large near the geometry walls and wake lines, and decrease gradually towards the far- 
field. In regions where these aspect ratios become smaller than unity, the streamwise point dis- 
tribution is coarsened by removing selected points, thus maintaining a nearly isotropic point 
distribution in the inviscid regions of flow. This filtered set of mesh points is then triangulated 
using a modified Delaunay triangulation criterion [10]. In its original form, a Delaunay tri- 
angulation of a given set of points tends to produce the most equiangular triangles possible, 
and is thus not well suited for the generation of highly stretched triangulations. The Delaunay 
criterion is thus modified according to the local stretching values. The local stretching value is 
used to define a stretched space, within which the local mesh-point distribution appears more 
isotropic. The Delaunay triangulation of the points in this stretched space is then constructed. 
This connected set of points is subsequently mapped back into physical space, thus producing 
the desired stretched triangulation. The actual construction of the stretched triangulation is per- 
formed in a two-stage process. A regular Delaunay triangulation of the mesh points in physi- 
cal space is initially constructed using Bowyer’s algorithm [14], This initial triangulation is 
conveniently used to smooth the distribution of stretchings throughout the flow-field by provid- 
ing the basis for discretizing a smoothing operator which is then applied to the stretching 
values. This Delaunay triangulation is then transformed into a modified (stretched) Delaunay 
triangulation, based on the smoothed local stretching distribution, by swapping the diagonals 
according to the Delaunay maximum minimum-angle criterion [10] applied in the stretched 
space. This criterion provides a basis for deciding whether or not an edge should be swapped 
by choosing the configuration which maximizes the smallest of the angles formed between the 
edge to be swapped and its immediate neighbors. 


t Supplied by L. Wigton, The Boeing Company. 
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A regular Delaunay triangulation of a set of points produces a smoothly varying mesh 
provided the mesh points are initially distributed evenly and isotropically. Similarly, a 
stretched or modified Delaunay triangulation of a set of points will result in a smoothly varying 
stretched mesh provided the mesh-point and the stretching distributions are closely coupled and 
provided the stretching distribution varies slowly with respect to the local mesh element size. 
The geometry-adaptive mesh-point distribution and the point filtering operation previously 
described, have thus been designed to ensure that these criteria are satisfied. 

Finally, the stretched triangulation can be smoothed in a postprocessing operation by 
slightly repositioning the points according to a smoothing operator discretized on the mesh. 
Once the points have been displaced, the mesh may no longer obey the modified Delaunay cri- 
terion, thus the edges may be swapped to recover this property. Multiple smoothing and 
edge-swapping passes may then be employed to further smooth the mesh. 

2.2. Adaptive Meshing Procedure 

Once the initial stretched unstructured mesh has been generated and the flow-field has 
been solved for on this mesh, a new adaptively refined mesh may be constructed by adding 
new points to the initial mesh in regions where large flow gradients or discretization errors are 
detected, and locally restructuring the mesh. In this work, the refinement criterion is based on 
the undivided difference of pressure and Mach number. Pressure gradients provide a good 
indication of inviscid flow phenomena, such as shocks and expansions, while Mach number 
variations can be used to identify viscous phenomena such as boundary layers and wakes. 
Although the mesh-point distribution can be refined in the adaptation process, the stretching 
distribution is held fixed. Thus, in order to maintain a close coupling between the final mesh- 
point distribution and the stretching distribution, an isotropic refinement strategy must be 
adopted. The variations of pressure and Mach number within each triangular element of the 
mesh are examined. When these are larger than some fraction of the average variations of 
pressure and Mach number over all cells of the mesh, three new mesh points are created, one 
at the midpoint of each of the three edges of the cell. Each new mesh point is assigned a 
stretching value taken as the average of the stretchings of the two points at either end of the 
generating mesh edge. The new points are then inserted into the existing mesh by locally res- 
tructuring the mesh according to Bowyer’s algorithm in the stretched space [14]. For each new 
mesh point, the triangles whose circumcircles are intersected by this new point are located. The 
union of all intersected triangles forms a convex polygonal region, which contains the new 
point. The existing structure in this region is removed, and a new structure is constructed by 
joining the new point to all the vertices of this polygonal region, as shown in Figure 2. By 
simultaneously refining all three sides of a mesh element, and by assigning average stretching 
values to the new points, directional biasing of the refinement process is avoided, and a smooth 
distribution of stretching is maintained. The use of Bowyer’s algorithm in the stretched space, 
which provides an effective method for constructing a stretched Delaunay triangulation through 
sequential point insertion and retriangulation, constitutes an ideal adaptive mesh enrichment 
strategy, and obviates the need for global mesh regeneration. 

When new boundary points are introduced, they must be repositioned onto the spline 
definition of the geometry boundary. For curved surfaces, this will not coincide with the mid- 
point of the original boundary mesh edge. For highly stretched meshes, the distance between 
these two locations may in fact be much larger than the average width of the elements in the 
vicinity of the boundary, and a restructuring of the entire region is required, as shown in Fig- 
ure 3. This is accomplished by drawing the line joining the mid-point of the boundary edge 
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being refined to the spline-displaced position of this new boundary point. The union of all tri- 
angular mesh elements intersected by this line, as well as elements whose circumcircles are 
intersected by the spline displaced boundary point, are tagged for reconstruction. This entire 
region is then restructured, as per the standard Bowyer algorithm, i.e. removing all mesh edges 
in this region, and forming a new structure by joining up the new mesh point to all the vertices 
of the polygonal region. 


3. FLOW SOLUTION 


In non-dimensional conservative vector form, the Navier-Stokes equations read 




+ V.F„ 


1 


V.F„ 


( 1 ) 


dt ‘ Re, 

where Re. denotes the overall flow Reynolds number, and w represents the conserved variables 


P 

P« 

w = pv (2) 

[ pE . 

P being the fluid density, u and v the cartesian velocity components, and E the internal energy. 
F c represents the convective flux vector, the components of which are algebraic functions of 
the conserved variables and the pressure, which itself can be related to the conserved variables 
through the perfect gas relation. F„ denotes the viscous flux vector, the components of which 
are functions of the first derivatives of the conserved variables. Equation (1) represents a set 
of partial differential equations which must be discretized in space in order to obtain a set of 
coupled ordinary differential equations, which can then be integrated in time to obtain the 
steady-state solution. Space discretization is performed using a Galerkin finite-element type for- 
mulation. Multiplying equation (1) by a test function <|>, and integrating over physical space 
yields 


dxdy + J|<t»V.F c dxdy = J|<t>V.F v dxdy (3) 

Integrating the flux integrals by parts, and neglecting boundary terms gives 

- - v * ‘My ( 4 ) 

In order to evaluate the flux balance equations at a vertex P, <|> is taken as a piecewise linear 
function which has the value 1 at node P, and vanishes at all other vertices. Therefore, the 
integrals in the above equation are non-zero only over triangles which contain the vertex P, 
thus defining the domain of influence of node P, as shown in Figure 4. To evaluate the above 
integrals, we make use of the fact that and are constant over a triangle, and evaluate spa- 
tial derivatives of <)> and w over a triangle using vertex values, by Green’s contour integral 
theorem. The convective fluxes F c are taken as piecewise linear functions in space, and the 
viscous fluxes F v are piecewise constant over each triangle, since they are formed from first 
derivatives in the flow variables. Evaluating the flux integrals with these assumptions, one 
obtains 


dt 


dxdy = ||F C ,V<|> dxdy 



dxdy 
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» F C A + F* 
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where the summations are over all triangles in the domain of influence, as shown in Figure 4. 
represents the directed (normal) edge length of the face of each triangle on the outer boun- 
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dary of the domain, F A F® are the convective fluxes at the two vertices at either end of this 
edge, and F* is the viscous flux in triangle e, e being a triangle in the domain of influence of 
<|>. If the integral on the left hand side of equation (5) is evaluated in the same manner, the time 
derivatives become coupled in space. Since we are not interested in the time-accuracy of the 
scheme, but only in the final steady-state solution, we employ the concept of a lumped mass 
matrix. This is equivalent to assuming w to be constant over the domain of influence while 
integrating the left hand side. Hence, we obtain 


Q. 


dw„ 


= z 


F A + F' 


AL 


dt % 2 *® Re, 

where the factor of 1/3 is introduced by the integration of <|> over the domain, and Q p 
represents the surface area of the domain of influence of P. For the convective fluxes, this 
procedure is equivalent to the vertex finite- volume formulation described in [4,5]. For a 
smoothly varying regular triangulation, the above formulation is second-order accurate. 


z i <f;-al«) 

«=i z 


( 6 ) 


Additional artificial dissipation terms are required to ensure stability and to capture 
shocks without producing numerical oscillations. This is necessary for both inviscid and 
viscous flow computations, since in the later case, large regions of the flow field behave essen- 
tially inviscidly and the physical viscosity is not sufficient to guarantee numerical stability for 
the type of mesh spacings typically employed. Artificial dissipation terms are thus constructed 
as a blend of a Laplacian and a biharmonic operator in the conserved flow variables. The 
Laplacian term represents a strong formally first-order accurate dissipation which is turned on 
only in the vicinity of a shock, and the biharmonic term represents a weaker second-order 
accurate dissipation which is employed in regions of smooth flow [5,11]. The spatially discre- 
tized equations are integrated in time to obtain the steady-state solution using a five-stage 
time-stepping scheme, where the convective terms are evaluated at each stage within a time 
step, and the dissipative terms (both physical and artificial) are only evaluated at the first, third, 
and fifth stages. This particular scheme has been designed to maintain stability in regions 
where the flow is dominated by viscous effects, and to rapidly dampen out high-frequency 
error components, which is an essential feature for a scheme intended to drive a multigrid algo- 
rithm. Convergence is accelerated by making use of local time-stepping, implicit residual 
averaging, and an unstructured multigrid algorithm [11]. 

The idea of a multigrid strategy is to accelerate the convergence to steady-state of a fine 
grid solution through corrections computed on coarser grids. An initial time step is performed 
on the fine grid, and the flow variables and residuals are then transferred to the coarse grid. A 
correction equation is constructed on the coarse grid by adding a forcing function to the origi- 
nal discretized equations. This forcing function is formed by taking the difference between the 
transferred residuals and the residuals of the transferred variables, thus ensuring that the evolu- 
tion of the coarse grid equations is driven by the fine grid residuals. Hence, when the fine grid 
residuals vanish, the coarse grid equations are identically satisfied, and generate zero correc- 
tions. After transferring values down from the fine grid, a time step is performed on the coarse 
grid, and the new values are transferred down to the next coarser grid. When the coarsest grid 
is reached, the computed corrections are successively interpolated back up to the finest grid, 
and the entire cycle is repeated. In the context of unstructured meshes, a sequence of coarse 
and fine meshes is best constructed by generating the individual meshes independently from 
one another (as opposed to subdividing a coarse mesh). Thus, in general, the coarse and fine 
meshes of a given sequence do not have any common mesh points or nested elements. Thus, 
the patterns for transferring the variables, residuals, and corrections back and forth between the 
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various meshes of the sequence must be determined in a preprocessing operation, where an 
efficient tree-search algorithm is employed [5]. 

Such a multigrid algorithm may be combined with an adaptive meshing strategy in a 
natural manner. First, a sequence of globally generated meshes is constructed, and multigrid 
time-stepping is performed on this sequence until a satisfactorily converged solution is 
obtained. At this point, a new adaptively refined mesh is generated, and the transfer patterns 
for transferring variables from the previous mesh to the new mesh are determined. The flow 
variables are then transferred to this new mesh, providing a starting solution, and multigrid 
time-stepping is resumed on this new sequence which now contains an additional fine mesh. 
The process may be repeated, as shown in Figure 5, each time adding a new finer mesh to the 
sequence, until a converged solution of the desired accuracy is obtained. 

4. TURBULENCE MODELING 
4.1. General Procedure 

The most widespread turbulence models in use currently are either of the multiple field- 
equation type, or of the algebraic type. While field-equation models (such as the K-e model) 
can be discretized and solved for on unstructured meshes in a straight-forward manner, the 
solution of additional field-equations can be considerably expensive, especially in near-wall 
regions where the equations may become very stiff numerically. Algebraic models, on the 
other hand, are relatively inexpensive to compute, and have demonstrated generally superior 
accuracy and reliability for limited classes of problems, such as high-Reynolds-number attached 
flows over streamlined bodies. However, such models typically require information concerning 
the distance of each mesh point from the nearest wall. Turbulence length scales, which are 
related to the local boundary-layer or wake thickness, are determined by scanning the appropri- 
ate flow values along specified streamwise stations. In the context of unstructured meshes, 
mesh points and thus flow variables do not naturally occur at regular streamwise locations. 
Thus, lines normal to the walls and viscous layers must be created, and flow variables interpo- 
lated onto these lines in order that turbulence length scales may be determined. This type of 
approach has previously been implemented for supersonic ramp geometries by Rostand [15]. 
However, in the present work, more complex geometries must be accommodated. A more 
sophisticated manner of generating normal mesh stations can be devised using structured 
hyperbolic mesh generation techniques [13]. Thus, a distribution of normal mesh lines which 
do not cross-over each other, and containing a smoothly varying normal distribution of points, 
can be obtained by generating a structured hyperbolic mesh about each geometry component, 
based on the boundary-point distribution of the global unstructured mesh on each component 
[12]. These normal mesh lines are terminated if they intersect a neighboring geometry com- 
ponent, thus ensuring that turbulence quantities in any given region of the flow-field are only 
dependent on the viscous layers and walls which are directly visible from that location (c.f. 
Figures 7 and 12). At each time-step in the flow solution phase, flow variables are interpolated 
onto the background turbulence mesh stations, and the Baldwin-Lomax [16] algebraic model is 
employed to compute eddy viscosity values along these stations. The standard Baldwin-Lomax 
model must be modified slightly, in order to enable the treatment of flows with multiple 
boundary-layers and wakes. Thus, the search for the maximum moment of vorticity, which is 
usually employed to determine a turbulence length scale, is limited to the region between the 
wall (or wake centerline), and the first point of vanishing vorticity off the wall (or wake center- 
line). Since a change in sign of the vorticity occurs in the region between two confluent shear 
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layers, turbulence length scales associated with vorticity fluctuations due to neighboring shear 
layers are ignored in this manner. The eddy viscosities are then interpolated back onto the 
unstructured mesh for subsequent use in the flow solver. In regions where multiple back- 
ground turbulence meshes overlap, the multiple eddy viscosity values (one from each mesh) 
interpolated back to the unstructured mesh are weighted by a factor proportional to the inverse 
of the distance from the corresponding wall, thus producing a smooth distribution of eddy 
viscosity throughout the flow-field. Flow variables can be repeatedly interpolated back and 
forth between the background turbulence mesh stations and the global unstructured mesh by 
computing and storing the interpolation weights and addresses in a preprocessing operation 
prior to the flow solution phase. This is accomplished by first triangulating the point distribu- 
tion of the local background turbulence meshes, and then making use of the same efficient 
search routines employed in the unstructured multigrid transfer process. 

4.2. Adaptive Meshing Procedure 

When an adaptive meshing strategy is employed, the background turbulence meshes must 
be adapted in a manner analogous to the refinement of the global unstructured mesh. A 
refinement field variable is thus constructed, which is assigned the value 1.0 in regions where 
the unstructured mesh is refined, and 0.0 in the remainder of the flow field. This variable is 
then interpolated onto the background turbulence meshes and used to determine the regions 
which require refinement. Since these background meshes have previously been triangulated 
for interpolation purposes, they may be refined by inserting new points and restructuring 
locally using Bowyer’s Delaunay triangulation algorithm [14J as described previously. How- 
ever, these new points must be added in a manner which preserves the original structure (that 
of normal mesh stations) of the background meshes. Hence, only new points along existing 
mesh stations are permitted, and when a new boundary point is inserted, an entire new tur- 
bulence mesh station extending up to the outer boundary of the background mesh is created 
[ 12 ]. 

After each adaptation process, the transfer patterns for interpolation between the newly 
refined global unstructured mesh and background meshes must be recomputed. In the context 
of the multigrid strategy, the turbulence model is only executed on the finest grid of the 
sequence. The computed eddy viscosity values are then interpolated up to the coarser unstruc- 
tured meshes where they are employed in the multigrid correction equations. The whole pro- 
cess is very efficient, and in general, the turbulence model is found to require less than 10% of 
the total time required within a single multigrid cycle. Memory requirements are, however, 
increased by close to 50% since extra variables and transfer coefficients must be stored for the 
background turbulence mesh stations. 

5. RESULTS 

5.1. Single Airfoil Geometry 

As an initial test case, the turbulent flow over an RAE 2822 airfoil has been computed. 
The freestream Mach number is 0.729, the Reynolds number is 6.5 million, the corrected 
incidence is 2.31 degrees, and transition is fixed at 0.03 chords. This constitutes a well docu- 
mented test case (Case 6) for turbulent transonic flow [17], which can be used to validate the 
present solver. The unstructured mesh employed for this case is depicted in Figure 6. This 
mesh contains 13,751 points of which 210 are on the airfoil surface. The average normal spac- 
ings of the triangles on the airfoil surface is 0.00001 chords, resulting in cell aspect ratios of 
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the order of 1000:1 near the wall. The background turbulence mesh stations employed for 
computing the algebraic turbulence model, which contain a total of 13,372 points, are depicted 
in Figure 7. The computed surface pressure distribution and skin friction distribution are 
displayed in Figures 8 and 9, respectively, where they are compared with experimental data 
from [17]. Both quantities are seen to compare favorably with the experimental results, and the 
computed lift coefficient of 0.7403 is well within the range reported in previously published 
computational solutions, using structured meshes [18]. A total of five meshes were employed 
in the multigrid sequence, with the coarsest mesh containing only 98 points. The convergence 
rate for this case, as measured by the decrease in the RMS average of the density residuals 
throughout the flow-field, versus the number of multigrid cycles, is depicted in Figure 10. An 
average residual reduction of 0.955 per multigrid cycle is achieved on the finest grid, resulting 
in a decrease of the residuals by 4 orders of magnitude over 200 cycles. Furthermore, the lift 
and drag coefficients were converged to four significant figures within 90 cycles. Since each 
multigrid cycle requires roughly 1.4 CPU seconds on a single processor of the CRAY-YMP 
computer, engineering solutions could thus be obtained in approximately 2 minutes for this 
case. 

5.2. Two-Element Airfoil Geometry 

The next test case consists of a main airfoil with a leading edge slat which has been the 
subject of extensive wind-tunnel tests, as part of a program aimed at determining the 
effectiveness of slats as maneuvering devices for fighter aircraft [19]. The test conditions con- 
sist of a freestream Mach number of 0.5, a chord Reynolds number of 4.5 million, and a 
corrected incidence of 7.5 degrees. At these conditions, the flow becomes supercritical, and a 
small shock is formed on the upper surface of the slat. The adapted mesh used to compute 
this flow is depicted in Figure 11. The mesh contains a total of 35,885 points, of which 418 
are on the surface of the main airfoil, and 421 on the surface of the slat. The minimum nor- 
mal spacing at the wall is 0.00001 chords, and cells of aspect ratios up to 1000:1 are observed. 
A total of 7 meshes were employed in the multigrid sequence with the last 3 meshes generated 
adaptively. Figure 12 illustrates the adaptively generated turbulence mesh stations employed by 
the algebraic turbulence model for this case. The computed Mach contours in the flow field 
are depicted in Figure 13 where a crisp resolution of the small localized shock provided by the 
adaptive meshing technique is observed. A good correlation between the computed and experi- 
mental surface pressure coefficients is displayed in Figure 14. This method predicts a small 
zone of separated flow behind the shock which has not been detected in previous calculations 
[19], while also providing a more accurate prediction of shock location and subsequent pres- 
sure recovery. The solution of this case required 15 minutes on a single processor of a 
CRAY-YMP, during which the fine grid residuals were reduced by 3 orders of magnitude over 
200 multigrid cycles, as shown in the convergence history plot of Figure 15. 

53 . Four-Element-Airfoil Geometry 

The final test case consists of a four-element airfoil configuration. This represents a truly 
complex geometry involving regions with sharp comers which is not easily amenable to struc- 
tured mesh techniques and is of considerable practical interest as it relates to the design of 
high-lift devices for commercial aircraft. A multigrid sequence of 6 meshes was employed to 
compute the flow over this configuration with the last 2 meshes generated adaptively. The 
finest mesh of the sequence, which contains a total of 48,691 points, is depicted in Figure 16. 
The mesh contains 243 points on the main airfoil, 327 points on the leading edge slat with 208 
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and 247 points located on the vane and flap respectively. The height of the smallest cells at 
the wall is of the order of 0.00002 chords for each airfoil element and cell aspect ratios up to 
500:1 are observed. The computed Mach contours for this case are depicted in Figure 17. The 
ffeestream Mach number is 0.1995, the chord Reynolds number is 1.187 million, and the 
corrected incidence is 16.02 degrees. At these conditions, the flow remains entirely subcritical. 
Compressibility effects are nevertheless important due to the large suction peaks generated 
about each airfoil. For example, in the suction peak on the upper surface of the leading-edge 
slat, the local Mach number achieves a value of 0.77. The computed surface pressure 
coefficients are compared with experimental wind tunnel data t in Figure 18, and good overall 
agreement, including the prediction of the height of the suction peaks, is observed. This case 
provides a good illustration of the importance of adaptive meshing in practical aerodynamic 
calculations. Adequate resolution of the strong suction peak on the upper surface of the slat can 
only be achieved with a very fine mesh resolution in this region. Failure to adequately capture 
this large suction peak results in the generation of numerical entropy which is then convected 
downstream, thus contaminating the solution in the downstream regions, and degenerating the 
global accuracy of the solution. Because these suction peaks are very localized, they are 
efficiently resolved with adaptive techniques. In order to obtain a similar resolution using glo- 
bal mesh refinement, of the order of 200,000 mesh points would be required, greatly increasing 
the cost of the computation. The convergence history for this case, as measured by the density 
residuals and the total lift coefficient versus the number of multigrid cycles, is depicted in Fig- 
ure 19. A total of 400 multigrid cycles were executed, which required roughly 35 minutes of 
single processor CRAY-YMP time, and 14 Mwords of memory. 

6. CONCLUSIONS 

A method for efficiently computing turbulent compressible flows over complex two- 
dimensional configurations has been presented. By employing fully unstructured meshes 
throughout the entire flow field, arbitrary geometries can be handled and adaptive meshing 
techniques may be extensively exploited. The method is efficient in that engineering solutions 
may generally be obtained in less than 100 multigrid cycles, and the turbulence model requires 
less than 10% of the total time required to compute a solution. Unstructured mesh algorithms 
incur substantial overheads due to the random nature of their data-sets and generally run up to 
three times slower than their structured counterparts on present day supercomputers due to the 
indirect addressing and gather- scatter operations which must be employed. However, adaptive 
meshing techniques enable the accurate resolution of complex flow-fields with a relatively 
small number of points easily outweighing the penalties incurred by the use of random data- 
sets. In the future, an adaptive meshing technique capable of modifying the local mesh stretch- 
ing, as well as the mesh-point distribution, should be pursued. Further turbulence modeling 
research is also required, either refining the current algebraic model or implementing an alter- 
nate field equation model for cases with large regions of separated flows. 
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Figure 1 

Overlapping Structured C-Meshes for Multi-Element Airfoil Geometry 
Providing Initial Definition of Mesh-Point Distribution and 
Local Stretching Factors 



a) 

Insertion of New Point 


Determination of Intersected Triangles 
and Removal of Mesh Structure in this Region 



c) 

Local Restructuring of Mesh in Affected Region 


Figure 2 

Illustration of Bowyer’s Algorithm for Delaunay Triangulation 
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Figure 4 

Domain of Influence of Finite-Element Basis Function and Equivalent 
Finite-Volume Control Volume 
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Figure 5 

Full Multigrid Algorithm Employed in Conjunction with Adaptive Meshing Strategy 



Figure 6 

Fully Unstructured Mesh with High Stretching Employed for Computing 
Turbulent Flow Past an RAE 2822 Airfoil (Number of Points = 13751) 
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Figure 7 

Turbulence Stations Employed for Computing Row Past an RAE 2822 Airfoil 
(Total Number of Points = 13372) 




Figure 8 

Comparison of Computed Surface Pressure with Experimental Measurements for Row 
over an RAE 2822 Airfoil (Mach = 0.729, Re=6.5 million. Incidence = 2.31 degrees) 
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Figure 9 

Comparison of Computed Skin Friction with Experimental Measurements for Flow 
over an RAE 2822 Airfoil (Mach — 0.729, Re=6.5 million. Incidence = 2.31 degrees) 
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Number of Cycles 

Nnode : 14045 Ncyct : 270 Ncycf : 200 

Resid 1 : 0.29E+01 Resid2 : 0.29E-03 Rate : 0.9548 


Figure 10 

Convergence Rate for Flow over an RAE 2822 Airfoil as Measured by the 
RMS average of the Density Residuals versus the Number of Multigrid Cycles 


No. of Supersonic Points 
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Figure 11 

Adaptively Generated Unstructured Mesh about Two-Element Airfoil 
Number of Nodes = 35,885 
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Figure 12 

Illustration of Turbulence Mesh Stations Employed in Algebraic Model 
Total Number of Points = 43,566 




